home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / lib / perl5 / Math / BigInt.pm < prev    next >
Text File  |  1995-07-02  |  10KB  |  348 lines
  1. package Math::BigInt;
  2.  
  3. %OVERLOAD = ( 
  4.                 # Anonymous subroutines:
  5. '+'    =>    sub {new BigInt &badd},
  6. '-'    =>    sub {new BigInt
  7.                $_[2]? bsub($_[1],${$_[0]}) : bsub(${$_[0]},$_[1])},
  8. '<=>'    =>    sub {new BigInt
  9.                $_[2]? bcmp($_[1],${$_[0]}) : bcmp(${$_[0]},$_[1])},
  10. 'cmp'    =>    sub {new BigInt
  11.                $_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
  12. '*'    =>    sub {new BigInt &bmul},
  13. '/'    =>    sub {new BigInt 
  14.                $_[2]? scalar bdiv($_[1],${$_[0]}) :
  15.              scalar bdiv(${$_[0]},$_[1])},
  16. '%'    =>    sub {new BigInt
  17.                $_[2]? bmod($_[1],${$_[0]}) : bmod(${$_[0]},$_[1])},
  18. '**'    =>    sub {new BigInt
  19.                $_[2]? bpow($_[1],${$_[0]}) : bpow(${$_[0]},$_[1])},
  20. 'neg'    =>    sub {new BigInt &bneg},
  21. 'abs'    =>    sub {new BigInt &babs},
  22.  
  23. qw(
  24. ""    stringify
  25. 0+    numify)            # Order of arguments unsignificant
  26. );
  27.  
  28. sub new {
  29.   my $foo = bnorm($_[1]);
  30.   die "Not a number initialized to BigInt" if $foo eq "NaN";
  31.   bless \$foo;
  32. }
  33. sub stringify { "${$_[0]}" }
  34. sub numify { 0 + "${$_[0]}" }    # Not needed, additional overhead
  35.                 # comparing to direct compilation based on
  36.                 # stringify
  37.  
  38. # arbitrary size integer math package
  39. #
  40. # by Mark Biggar
  41. #
  42. # Canonical Big integer value are strings of the form
  43. #       /^[+-]\d+$/ with leading zeros suppressed
  44. # Input values to these routines may be strings of the form
  45. #       /^\s*[+-]?[\d\s]+$/.
  46. # Examples:
  47. #   '+0'                            canonical zero value
  48. #   '   -123 123 123'               canonical value '-123123123'
  49. #   '1 23 456 7890'                 canonical value '+1234567890'
  50. # Output values always always in canonical form
  51. #
  52. # Actual math is done in an internal format consisting of an array
  53. #   whose first element is the sign (/^[+-]$/) and whose remaining 
  54. #   elements are base 100000 digits with the least significant digit first.
  55. # The string 'NaN' is used to represent the result when input arguments 
  56. #   are not numbers, as well as the result of dividing by zero
  57. #
  58. # routines provided are:
  59. #
  60. #   bneg(BINT) return BINT              negation
  61. #   babs(BINT) return BINT              absolute value
  62. #   bcmp(BINT,BINT) return CODE         compare numbers (undef,<0,=0,>0)
  63. #   badd(BINT,BINT) return BINT         addition
  64. #   bsub(BINT,BINT) return BINT         subtraction
  65. #   bmul(BINT,BINT) return BINT         multiplication
  66. #   bdiv(BINT,BINT) return (BINT,BINT)  division (quo,rem) just quo if scalar
  67. #   bmod(BINT,BINT) return BINT         modulus
  68. #   bgcd(BINT,BINT) return BINT         greatest common divisor
  69. #   bnorm(BINT) return BINT             normalization
  70. #
  71.  
  72. $zero = 0;
  73.  
  74.  
  75. # normalize string form of number.   Strip leading zeros.  Strip any
  76. #   white space and add a sign, if missing.
  77. # Strings that are not numbers result the value 'NaN'.
  78.  
  79. sub bnorm { #(num_str) return num_str
  80.     local($_) = @_;
  81.     s/\s+//g;                           # strip white space
  82.     if (s/^([+-]?)0*(\d+)$/$1$2/) {     # test if number
  83.     substr($_,$[,0) = '+' unless $1; # Add missing sign
  84.     s/^-0/+0/;
  85.     $_;
  86.     } else {
  87.     'NaN';
  88.     }
  89. }
  90.  
  91. # Convert a number from string format to internal base 100000 format.
  92. #   Assumes normalized value as input.
  93. sub internal { #(num_str) return int_num_array
  94.     local($d) = @_;
  95.     ($is,$il) = (substr($d,$[,1),length($d)-2);
  96.     substr($d,$[,1) = '';
  97.     ($is, reverse(unpack("a" . ($il%5+1) . ("a5" x ($il/5)), $d)));
  98. }
  99.  
  100. # Convert a number from internal base 100000 format to string format.
  101. #   This routine scribbles all over input array.
  102. sub external { #(int_num_array) return num_str
  103.     $es = shift;
  104.     grep($_ > 9999 || ($_ = substr('0000'.$_,-5)), @_);   # zero pad
  105.     &bnorm(join('', $es, reverse(@_)));    # reverse concat and normalize
  106. }
  107.  
  108. # Negate input value.
  109. sub bneg { #(num_str) return num_str
  110.     local($_) = &bnorm(@_);
  111.     vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0';
  112.     s/^H/N/;
  113.     $_;
  114. }
  115.  
  116. # Returns the absolute value of the input.
  117. sub babs { #(num_str) return num_str
  118.     &abs(&bnorm(@_));
  119. }
  120.  
  121. sub abs { # post-normalized abs for internal use
  122.     local($_) = @_;
  123.     s/^-/+/;
  124.     $_;
  125. }
  126.  
  127. # Compares 2 values.  Returns one of undef, <0, =0, >0. (suitable for sort)
  128. sub bcmp { #(num_str, num_str) return cond_code
  129.     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  130.     if ($x eq 'NaN') {
  131.     undef;
  132.     } elsif ($y eq 'NaN') {
  133.     undef;
  134.     } else {
  135.     &cmp($x,$y);
  136.     }
  137. }
  138.  
  139. sub cmp { # post-normalized compare for internal use
  140.     local($cx, $cy) = @_;
  141.     $cx cmp $cy
  142.     &&
  143.     (
  144.     ord($cy) <=> ord($cx)
  145.     ||
  146.     ($cx cmp ',') * (length($cy) <=> length($cx) || $cy cmp $cx)
  147.     );
  148. }
  149.  
  150. sub badd { #(num_str, num_str) return num_str
  151.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  152.     if ($x eq 'NaN') {
  153.     'NaN';
  154.     } elsif ($y eq 'NaN') {
  155.     'NaN';
  156.     } else {
  157.     @x = &internal($x);             # convert to internal form
  158.     @y = &internal($y);
  159.     local($sx, $sy) = (shift @x, shift @y); # get signs
  160.     if ($sx eq $sy) {
  161.         &external($sx, &add(*x, *y)); # if same sign add
  162.     } else {
  163.         ($x, $y) = (&abs($x),&abs($y)); # make abs
  164.         if (&cmp($y,$x) > 0) {
  165.         &external($sy, &sub(*y, *x));
  166.         } else {
  167.         &external($sx, &sub(*x, *y));
  168.         }
  169.     }
  170.     }
  171. }
  172.  
  173. sub bsub { #(num_str, num_str) return num_str
  174.     &badd($_[$[],&bneg($_[$[+1]));    
  175. }
  176.  
  177. # GCD -- Euclids algorithm Knuth Vol 2 pg 296
  178. sub bgcd { #(num_str, num_str) return num_str
  179.     local($x,$y) = (&bnorm($_[$[]),&bnorm($_[$[+1]));
  180.     if ($x eq 'NaN' || $y eq 'NaN') {
  181.     'NaN';
  182.     } else {
  183.     ($x, $y) = ($y,&bmod($x,$y)) while $y ne '+0';
  184.     $x;
  185.     }
  186. }
  187.  
  188. # routine to add two base 1e5 numbers
  189. #   stolen from Knuth Vol 2 Algorithm A pg 231
  190. #   there are separate routines to add and sub as per Kunth pg 233
  191. sub add { #(int_num_array, int_num_array) return int_num_array
  192.     local(*x, *y) = @_;
  193.     $car = 0;
  194.     for $x (@x) {
  195.     last unless @y || $car;
  196.     $x -= 1e5 if $car = (($x += shift(@y) + $car) >= 1e5);
  197.     }
  198.     for $y (@y) {
  199.     last unless $car;
  200.     $y -= 1e5 if $car = (($y += $car) >= 1e5);
  201.     }
  202.     (@x, @y, $car);
  203. }
  204.  
  205. # subtract base 1e5 numbers -- stolen from Knuth Vol 2 pg 232, $x > $y
  206. sub sub { #(int_num_array, int_num_array) return int_num_array
  207.     local(*sx, *sy) = @_;
  208.     $bar = 0;
  209.     for $sx (@sx) {
  210.     last unless @y || $bar;
  211.     $sx += 1e5 if $bar = (($sx -= shift(@sy) + $bar) < 0);
  212.     }
  213.     @sx;
  214. }
  215.  
  216. # multiply two numbers -- stolen from Knuth Vol 2 pg 233
  217. sub bmul { #(num_str, num_str) return num_str
  218.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  219.     if ($x eq 'NaN') {
  220.     'NaN';
  221.     } elsif ($y eq 'NaN') {
  222.     'NaN';
  223.     } else {
  224.     @x = &internal($x);
  225.     @y = &internal($y);
  226.     &external(&mul(*x,*y));
  227.     }
  228. }
  229.  
  230. # multiply two numbers in internal representation
  231. # destroys the arguments, supposes that two arguments are different
  232. sub mul { #(*int_num_array, *int_num_array) return int_num_array
  233.     local(*x, *y) = (shift, shift);
  234.     local($signr) = (shift @x ne shift @y) ? '-' : '+';
  235.     @prod = ();
  236.     for $x (@x) {
  237.       ($car, $cty) = (0, $[);
  238.       for $y (@y) {
  239.     $prod = $x * $y + $prod[$cty] + $car;
  240.     $prod[$cty++] =
  241.       $prod - ($car = int($prod * 1e-5)) * 1e5;
  242.       }
  243.       $prod[$cty] += $car if $car;
  244.       $x = shift @prod;
  245.     }
  246.     ($signr, @x, @prod);
  247. }
  248.  
  249. # modulus
  250. sub bmod { #(num_str, num_str) return num_str
  251.     (&bdiv(@_))[$[+1];
  252. }
  253.  
  254. sub bdiv { #(dividend: num_str, divisor: num_str) return num_str
  255.     local (*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  256.     return wantarray ? ('NaN','NaN') : 'NaN'
  257.     if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0');
  258.     return wantarray ? ('+0',$x) : '+0' if (&cmp(&abs($x),&abs($y)) < 0);
  259.     @x = &internal($x); @y = &internal($y);
  260.     $srem = $y[$[];
  261.     $sr = (shift @x ne shift @y) ? '-' : '+';
  262.     $car = $bar = $prd = 0;
  263.     if (($dd = int(1e5/($y[$#y]+1))) != 1) {
  264.     for $x (@x) {
  265.         $x = $x * $dd + $car;
  266.         $x -= ($car = int($x * 1e-5)) * 1e5;
  267.     }
  268.     push(@x, $car); $car = 0;
  269.     for $y (@y) {
  270.         $y = $y * $dd + $car;
  271.         $y -= ($car = int($y * 1e-5)) * 1e5;
  272.     }
  273.     }
  274.     else {
  275.     push(@x, 0);
  276.     }
  277.     @q = (); ($v2,$v1) = @y[-2,-1];
  278.     while ($#x > $#y) {
  279.     ($u2,$u1,$u0) = @x[-3..-1];
  280.     $q = (($u0 == $v1) ? 99999 : int(($u0*1e5+$u1)/$v1));
  281.     --$q while ($v2*$q > ($u0*1e5+$u1-$q*$v1)*1e5+$u2);
  282.     if ($q) {
  283.         ($car, $bar) = (0,0);
  284.         for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
  285.         $prd = $q * $y[$y] + $car;
  286.         $prd -= ($car = int($prd * 1e-5)) * 1e5;
  287.         $x[$x] += 1e5 if ($bar = (($x[$x] -= $prd + $bar) < 0));
  288.         }
  289.         if ($x[$#x] < $car + $bar) {
  290.         $car = 0; --$q;
  291.         for ($y = $[, $x = $#x-$#y+$[-1; $y <= $#y; ++$y,++$x) {
  292.             $x[$x] -= 1e5
  293.             if ($car = (($x[$x] += $y[$y] + $car) > 1e5));
  294.         }
  295.         }   
  296.     }
  297.     pop(@x); unshift(@q, $q);
  298.     }
  299.     if (wantarray) {
  300.     @d = ();
  301.     if ($dd != 1) {
  302.         $car = 0;
  303.         for $x (reverse @x) {
  304.         $prd = $car * 1e5 + $x;
  305.         $car = $prd - ($tmp = int($prd / $dd)) * $dd;
  306.         unshift(@d, $tmp);
  307.         }
  308.     }
  309.     else {
  310.         @d = @x;
  311.     }
  312.     (&external($sr, @q), &external($srem, @d, $zero));
  313.     } else {
  314.     &external($sr, @q);
  315.     }
  316. }
  317.  
  318. # compute power of two numbers -- stolen from Knuth Vol 2 pg 233
  319. sub bpow { #(num_str, num_str) return num_str
  320.     local(*x, *y); ($x, $y) = (&bnorm($_[$[]), &bnorm($_[$[+1]));
  321.     if ($x eq 'NaN') {
  322.     'NaN';
  323.     } elsif ($y eq 'NaN') {
  324.     'NaN';
  325.     } elsif ($x eq '+1') {
  326.     '+1';
  327.     } elsif ($x eq '-1') {
  328.     &bmod($x,2) ? '-1': '+1';
  329.     } elsif ($y =~ /^-/) {
  330.     'NaN';
  331.     } elsif ($x eq '+0' && $y eq '+0') {
  332.     'NaN';
  333.     } else {
  334.     @x = &internal($x);
  335.     local(@pow2)=@x;
  336.     local(@pow)=&internal("+1");
  337.     local($y1,$res,@tmp1,@tmp2)=(1); # need tmp to send to mul
  338.     while ($y ne '+0') {
  339.       ($y,$res)=&bdiv($y,2);
  340.       if ($res ne '+0') {@tmp=@pow2; @pow=&mul(*pow,*tmp);}
  341.       if ($y ne '+0') {@tmp=@pow2;@pow2=&mul(*pow2,*tmp);}
  342.     }
  343.     &external(@pow);
  344.     }
  345. }
  346.  
  347. 1;
  348.